1 Library Settings

Note: we must prevent the scientific notation issue. For example, when using paste("abc", 10000), we want the results to be "abc 10000" instead of "abc 1e+05". Hence, using options("scipen" =10) would deal with this issue. The reason of doing this is that the solve_POMDP() would first call write_POMDP() which convert the transition probabilities, observation probabilities and rewards to be a file.POMDP. The key of method is using paste() function. If there is no prevention for scientific notation, then the results from paste() would be some characters containing "1e-09", thereby causing that the Initilization POMDP part have several syntax error.

2 Introduction

In this R markdown, we will redo the development of cyber vulnerabilities policies based on POMDP, which is in our submitted ASMBI paper. The parameters of POMDP have been estimated and shown in the paper.

3 POMDP model

3.2 Transitions Prob., Observations Prob. and Rewards

3.2.2 Data structure for POMDP solver

Next, we modify the data structure for POMDP solver.

(1) Transition Probabilities

transition_prob is a \(OS\times A\times S\times S'\) tensor.

## $`auto-patching`
##        [,1]   [,2]   [,3]   [,4]
## [1,] 0.9171 0.0800 0.0020 0.0009
## [2,] 0.0031 0.9905 0.0059 0.0005
## [3,] 0.0000 0.0032 0.9905 0.0063
## [4,] 0.0000 0.0000 0.0032 0.9968
## 
## $`research-accept`
##        [,1]   [,2]   [,3]   [,4]
## [1,] 0.9963 0.0037 0.0000 0.0000
## [2,] 0.4299 0.5664 0.0037 0.0000
## [3,] 0.0316 0.3983 0.5664 0.0037
## [4,] 0.0091 0.0604 0.0205 0.9100
## 
## $`research-compensate`
##        [,1]   [,2]   [,3] [,4]
## [1,] 1.0000 0.0000 0.0000    0
## [2,] 1.0000 0.0000 0.0000    0
## [3,] 0.0315 0.9648 0.0037    0
## [4,] 0.0091 0.0604 0.9305    0
## 
## $remediation
##      [,1] [,2] [,3] [,4]
## [1,]  0.5  0.5    0    0
## [2,]  0.5  0.5    0    0
## [3,]  0.5  0.5    0    0
## [4,]  0.5  0.5    0    0

(2) Observation Probabilities

observation_prob is a \(OS\times Criticality\times A\times S'\times O\).

## $`auto-patching`
##           [,1]      [,2]      [,3]      [,4]      [,5]
## [1,] 0.3582727 0.6417273 0.0000000 0.0000000 0.0000000
## [2,] 0.3582727 0.0000000 0.6417273 0.0000000 0.0000000
## [3,] 0.3582727 0.0000000 0.0000000 0.6417273 0.0000000
## [4,] 0.3582727 0.0000000 0.0000000 0.0000000 0.6417273
## 
## $`research-accept`
##           [,1]      [,2]      [,3]      [,4]      [,5]
## [1,] 0.3582727 0.6417273 0.0000000 0.0000000 0.0000000
## [2,] 0.3582727 0.0000000 0.6417273 0.0000000 0.0000000
## [3,] 0.3582727 0.0000000 0.0000000 0.6417273 0.0000000
## [4,] 0.3582727 0.0000000 0.0000000 0.0000000 0.6417273
## 
## $`research-compensate`
##           [,1]      [,2]      [,3]      [,4]      [,5]
## [1,] 0.3582727 0.6417273 0.0000000 0.0000000 0.0000000
## [2,] 0.3582727 0.0000000 0.6417273 0.0000000 0.0000000
## [3,] 0.3582727 0.0000000 0.0000000 0.6417273 0.0000000
## [4,] 0.3582727 0.0000000 0.0000000 0.0000000 0.6417273
## 
## $remediation
##           [,1]      [,2]      [,3]      [,4]      [,5]
## [1,] 0.3582727 0.6417273 0.0000000 0.0000000 0.0000000
## [2,] 0.3582727 0.0000000 0.6417273 0.0000000 0.0000000
## [3,] 0.3582727 0.0000000 0.0000000 0.6417273 0.0000000
## [4,] 0.3582727 0.0000000 0.0000000 0.0000000 0.6417273

Here, we consider the error scenarios to vary observation matrix. For example, let error be 10%.

  • For actual \(s' = Low\), observing medium, high or critical state is impossible. So such error can not affect this.
  • For actual \(s' = Medium\), we might observe \(Low\) with 10% probability. Then,

\[ P_{error}(o = Medium|s' = Medium) = P(o = Medium|s'= Medium)*(1-10\%) \\ P_{error}(o = Low|s' = Medium) = P(o = Medium|s'= Medium)*(10\%) \]

  • For actual \(s' = High\), we might observe \(Low\) or \(Medium\) with 10% probability. Then,

\[ P_{error}(o = High|s' = High) = P(o = High|s'= High)*(1-10\%) \\ P_{error}(o = Medium|s' = High) = P(o = High|s'= High)*(10\%/2)\\ P_{error}(o = Low|s' = High) = P(o = High|s'= High)*(10\%/2) \]

  • For actual \(s' = High\), we might observe \(Low\) or \(Medium\) with 10% probability. Then,

\[ P_{error}(o = High|s' = High) = P(o = High|s'= High)*(1-10\%) \\ P_{error}(o = Medium|s' = High) = P(o = High|s'= High)*(10\%/2)\\ P_{error}(o = Low|s' = High) = P(o = High|s'= High)*(10\%/2) \]

  • For actual \(s' = Critical\), we might observe \(Low\), \(Medium\) or \(Critical\) with 10% probability. Then,

\[ P_{error}(o = Critical|s' = Critical) = P(o = Critical|s'= Critical)*(1-10\%) \\ P_{error}(o = High|s' = Critical) = P(o = Critical|s'= Critical)*(10\%/3) \\ P_{error}(o = Medium|s' = Critical) = P(o = Critical|s'= Critical)*(10\%/3)\\ P_{error}(o = Low|s' = Critical) = P(o = Critical|s'= Critical)*(10\%/3) \]

The above procedure has been put into the following function Consider_Error_Observation

(3) Reward

reward is a \(OS\times Criticality\times Admin \times Dataframe\) where \(Dataframe\) contains columns: action, start-state, end-state, observation and reward.

4 Solving POMDPs Functions

There are some issues in solve_POMDP() from pomdp package. Specifically, the issues happen in $pg. Therefore, we have to modify the solve_POMDP() function. In addition, some functions have been made.

4.2 belief_update function

Every time agent taking an action \(a\) and observes \(o\), its belief is updated by Bayes rule:

\[ b^{ao}(s') = \frac{p(o|s',a)}{p(o|b,a)}\sum\limits_{s\in S}p(s'|s,a)*b(s) \]

where \(p(s'|s,a)\) and \(p(o|s',a)\) are defined by the model parameters transition probability and observation probability, and

\[ p(o|b,a) = \sum\limits_{s'\in S}p(o|s', a)\sum\limits_{s\in S}p(s'|s,a)*b(s) \]

whichh is normalization term.

Here, we define belief_update() function to implement the above formulation.

belief_update <- function(b_0, model, act_index, obs_index){
  if (class(model)!= "POMDP_model"){
    stop("the 'model' argument must be constructed by POMDP() funtion")
  }else {
    n_observations <- length(model$observations) # number of observations
    n_state <- length(model$states) # number of states
    # if (act_index > length(model$actions)){
    #   stop("Error in act_index: it is out of bound")
    # }
    # if(obs_index > n_observations){
    #   stop("Error in obs_index: it is out of bound")
    # }
       transition_prob_act <- model$transition_prob[[act_index]]%>%as.matrix() # transition prob of action
       observation_prob_act <- model$observation_prob[[act_index]]%>%as.matrix() #observation prob of action
       if (observation_prob_act[1,1] == "uniform"){
         observation_prob_act <- matrix(1/n_observations, nrow = n_state, ncol = n_observations)
       }
       # compute the denominator (normalization term)
       sum_1<-0
       sum_2<-0
       for (s_end in 1:n_state){
         for (s_start in 1:n_state){
           x<-b_0[s_start]*transition_prob_act[s_start, s_end]
           sum_1<-sum_1+x
           }
         y<-sum_1*observation_prob_act[s_end, obs_index]
         sum_2<-sum_2+y
         #print(sum_1)#debug
         }
       p_normalizer<-sum_2
       # update belief
       if(p_normalizer == 0){
           b_update <- rep(NA_real_, n_state)
           cat(paste("Warnings in ", model$name,": The belief update is NAs because it is impossible to observe", 
                     model$observations[obs_index], 
                     "observation under action", 
                     model$actions[act_index],
                     "when the start belief is (",
                     paste(b_0, collapse = ","), ")",
                     sep = " "))
           return(b_update)
           }else{
             b_update<-b_0
             sum<-0
             for (s_end in 1:n_state){
               for (s_start in 1:n_state){
                 x<-b_0[s_start]*transition_prob_act[s_start, s_end]
                 sum<-sum+x
                 }
               b_update[s_end]<-observation_prob_act[s_end, obs_index]*sum/p_normalizer
               }
             return(b_update)
       }
    }
}

#belief_update(b_0 = c(1,0,0,0), model = model_POMDPs[[1+12]], act_index = 1, obs_index = 3)

4.3 create_PG function

Since the pg dataframe has some problem, here we are using create_PG function to re-construct the policy graph.

The input of this function is:

  • model: it is POMDP model which has been constructed by POMDP() function

  • alpha: it is a dataframe that contains alpha vectors. Each row refers to an alpha vector. The number of columns is equal to the number of states.

  • action_index_alpha: it is a vector of actions index. Each action is associated with the row of alpha dataframe.

  • belief_matrix: it is a matrix. Each row refers to a vector of belief, for example, (0.25, 0.25, 0.25, 0.25). The number of columns is equal to the number of states.

create_PG <- function(model, alpha, action_index_alpha, belief_matrix){
  if (class(model)!= "POMDP_model"){
    stop("Error in model: the 'model' argument must be constructed by POMDP() funtion")
  }else {
    if(is.null(belief_matrix)){
      pg <- NULL
      return(pg)
    }else{
      alpha_matrix <- as.matrix(alpha) 
      n_alpha <-dim(alpha_matrix)[1]
      if(length(action_index_alpha)!=n_alpha){
        stop("Error in action_index_alpha: the length of this vector has to be the same as the number of rows of alpha")
      }
      # PG construction: 
      # initialize optimal action index for each row of belief matrix
      act_opt_index <- rep(0, dim(belief_matrix)[1])
      # initialize PG
      n_observations <- length(model$observations) 
      PG_int <- matrix(0, nrow = dim(belief_matrix)[1], ncol = 2 + n_observations)
      
      # for each row of belief matrix, find the index of alpha vector that gives maximum value 
      # and then find the opt action index
      for (i in 1:dim(belief_matrix)[1]){
        if (is.na(belief_matrix[i, 1])){ #if the belief points are NAs
          PG_int[i, 1] <-NA_integer_
        }else{
          alpha_opt_index <- which.max(alpha_matrix %*% belief_matrix[i,])
          act_opt_index[i] <- action_index_alpha[alpha_opt_index]
          PG_int[i,1] <- alpha_opt_index
        }
        
      }
      
      PG_int[ ,2] <- act_opt_index 

      # for each row of belief matrix, get the alpha opt index for each observation
      for (i in 1:dim(belief_matrix)[1]){
        for (j in 1:n_observations){
          if (is.na(belief_matrix[i, 1])){ #if the belief points are NAs
            PG_int[i, j+2] <- NA_real_
          }else{
            b_update <- belief_update(b_0 = belief_matrix[i, ],
                                    model = model,
                                    act_index = act_opt_index[i],
                                    obs_index = j)
            if (is.na(b_update[1])){ #if the updated belief is NA
              PG_int[i, j+2] <- NA_real_
            }else{
              PG_int[i, j+2] <- which.max(alpha_matrix %*% b_update)
              }
            }
          }
        }
      #rename and construct PG as dataframe
      pg <- as.data.frame(PG_int)
      # renamig the columns
      colnames(pg)[1] <- "belief_state"
      colnames(pg)[2] <- "action"
      for (i in 1:n_observations) {
        colnames(pg)[i+2] <- model$observations[i]
        }
      # renaming the actions
      for (i in 1:length(model$actions)) {
        pg[pg[,2]==i,2] <- model$actions[i]
      }

      return(pg)
    }
  }
}

# create_PG(model = LinuxProblem2,
#           alpha = Linux_sol$alpha,
#           action_index_alpha = action_index_alpha,
#           belief_matrix = belief_prop_matrix)

4.4 create_PG2 function

This one will generate PG with the action name or action index

create_PG2 <- function(model, alpha, action_index_alpha, belief_matrix,
                       get_action_name = TRUE){
  if (class(model)!= "POMDP_model"){
    stop("Error in model: the 'model' argument must be constructed by POMDP() funtion")
  }else {
    if(is.null(belief_matrix)){
      pg <- NULL
      return(pg)
    }else{
      alpha_matrix <- as.matrix(alpha) 
      n_alpha <-dim(alpha_matrix)[1]
      if(length(action_index_alpha)!=n_alpha){
        stop("Error in action_index_alpha: the length of this vector has to be the same as the number of rows of alpha")
      }
      # PG construction: 
      # initialize optimal action index for each row of belief matrix
      act_opt_index <- rep(0, dim(belief_matrix)[1])
      # initialize PG
      n_observations <- length(model$observations) 
      PG_int <- matrix(0, nrow = dim(belief_matrix)[1], ncol = 2 + n_observations)
      
      # for each row of belief matrix, find the index of alpha vector that gives maximum value 
      # and then find the opt action index
      for (i in 1:dim(belief_matrix)[1]){
        if (is.na(belief_matrix[i, 1])){ #if the belief points are NAs
          PG_int[i, 1] <-NA_integer_
        }else{
          alpha_opt_index <- which.max(alpha_matrix %*% belief_matrix[i,])
          act_opt_index[i] <- action_index_alpha[alpha_opt_index]
          PG_int[i,1] <- i-1
        }
        
      }
      
      if (get_action_name){
        PG_int[ ,2] <- model$actions[act_opt_index] 
      }else{
        PG_int[ ,2] <- act_opt_index
      }
      

      # for each row of belief matrix, get the alpha opt index for each observation
      for (i in 1:dim(belief_matrix)[1]){
        for (j in 1:n_observations){
          if (is.na(belief_matrix[i, 1])){ #if the belief points are NAs
            if(get_action_name){
              PG_int[i, j+2] <- NA_character_
            }else{
               PG_int[i, j+2] <- NA_integer_
            }
            
          }else{
            b_update <- belief_update(b_0 = belief_matrix[i, ],
                                    model = model,
                                    act_index = act_opt_index[i],
                                    obs_index = j)
            if (is.na(b_update[1])){ #if the updated belief is NA
              if(get_action_name){
                PG_int[i, j+2] <- NA_character_
              }else{
                PG_int[i, j+2] <- NA_integer_
               }
            }else{
              alpha_opt_index <- which.max(alpha_matrix %*% b_update)
              act_index <- action_index_alpha[alpha_opt_index] # get action index
              if (get_action_name){
                PG_int[i, j+2] <- model$actions[act_index]
              }else{
                PG_int[i, j+2] <- act_index
              }
              }
            }
          }
        }
      #rename and construct PG as dataframe
      pg <- as.data.frame(PG_int)
      # renamig the columns
      colnames(pg)[1] <- "belief_scenario"
      colnames(pg)[2] <- "action"
      for (i in 1:n_observations) {
        colnames(pg)[i+2] <- model$observations[i]
        }
      # # renaming the actions
      # for (i in 1:length(model$actions)) {
      #   pg[pg[,2]==paste(i),2] <- model$actions[i]
      # }

      return(pg)
    }
  }
}

# create_PG2(model = model_POMDPs_solved[[5]]$model,
#           alpha = model_POMDPs_solved[[5]]$solution$alpha,
#           action_index_alpha = model_POMDPs_solved[[5]]$solution$action_index_alpha,
#            belief_matrix = matrix(c(0.25, 0.25, 0.25, 0.25,
#                                     1, 0, 0, 0,
#                                     0, 1, 0, 0,
#                                     0, 0, 1, 0,
#                                     0, 0, 0, 1), nrow = 5, ncol = 4, byrow = TRUE),
#           get_action_name = TRUE)

4.5 solve_POMDP2 function

This function is modified based on solve_POMDP. The key difference is that we re-calculate the pg dataframe based on belief_update function, create_PG function, and alpha vector instead of reading the .pg file from the results of pomdp-solve.exe.

solve_POMDP2 <- function(
  model,
  horizon = NULL,
  method = "grid",
  parameter = NULL,
  verbose = FALSE) {

  ### DEBUG
  #verbose <- TRUE
  ###
    
  methods <- c("grid", "enum", "twopass", "witness", "incprune")
  ### Not implemented:  "linsup", "mcgs"
  method <- match.arg(method, methods)
  
  ### write model to file
  file_prefix <- tempfile(pattern = "model")
  pomdp_filename <- paste0(file_prefix, ".POMDP") 
  pomdp::write_POMDP(model, pomdp_filename)
    
  ### running the POMDP code
  
  if(!is.null(parameter)) {
    paras <- sapply(names(parameter), FUN = function(n) paste0("-", n, " ", parameter[[n]]))
  } else paras <- ""
  
  solver_output <- system2(find_pomdpsolve(), #invoke a system command
    args = c(paste("-pomdp", pomdp_filename),
             paste("-method", method),
             (if(!is.null(horizon)) paste("-horizon", horizon) else ""),
             paras, 
             "-fg_save true"),
    stdout = TRUE, 
    stderr = TRUE, 
    wait = TRUE
  )
    
  ## the verbose mode: printing all the outputs from pomdp solver
  if(verbose) cat(paste(solver_output, "\n"))
  
  ### importing the outputs and results (pomdp-solve adds the PID to the file prefix)
  file_prefix <- gsub("^o = (.*)\\s*$", "\\1", solver_output[16]) 
  
  ## Creating result files' names and extensions
  pomdp_filename <- paste0(file_prefix, ".POMDP") 
  pg_filename <- paste0(file_prefix, ".pg")
  belief_filename <- paste0(file_prefix, ".belief")
  alpha_filename <- paste0(file_prefix, ".alpha")
  
  ## importing alpha file
  number_of_states <- length(model$states)
  number_of_observations <- length(model$observations)
  number_of_actions <- length(model$actions)
  observations <- model$observations
  actions <- model$actions
  states <- model$states
  start <- model$start
  
  alpha <- read.table(alpha_filename, header = FALSE, sep = "\n")
  alpha <- as.matrix(alpha)
  toDel <- seq(1,dim(alpha)[1],2) #the row index which contains the value of action
  action_index_alpha <- alpha[toDel, 1]%>%as.numeric() # the actions index associated with each alpha vector
  action_index_alpha <- action_index_alpha +1
  alpha <- alpha[-toDel,]
  alpha <- unlist(strsplit(alpha, " "))
  alpha<- matrix(as.numeric(alpha), ncol = number_of_states, byrow = TRUE)
  alpha_matrix <- alpha
  alpha <- as.data.frame(alpha)
  
  # renaming the columns
  for (i in 1:number_of_states) {
    colnames(alpha)[i] <- paste("coeffecient", i)
  }
  
  ## importing pg file
  # pg <- read.table(pg_filename, header = FALSE, sep = " ", quote = "\"", 
  #   dec = ".", na.strings = NA, numerals = "no.loss")
  # pg <- as.matrix(pg)
  # pg <- pg[,c(-3,-dim(pg)[2])] #there 2 NA columns that need to be removed
  # pg <- pg+1 #index has to start from 1 not 0
  # pg_matrix <- pg
  # pg <- as.data.frame(pg)
  # if (dim(pg)[2]==1 ) {
  #   pg <- t(pg)
  #   pg <-as.data.frame(pg)
  # }
  
  
  ## importing belief file if it exists
  if(file.exists(belief_filename)) {
    belief <- read.table(belief_filename) 
    belief_matrix <- as.matrix(belief)
    belief <- as.data.frame(belief_matrix)
    
    # renaming the columns
    for (i in 1:number_of_states) {
      colnames(belief)[i] <- states[i]
    }
  
    ## finding the respective proportions for each line (node)
    for (i in 1:dim(belief_matrix)[1]) {
      belief[i,number_of_states+1]<- which.max(alpha_matrix %*% belief_matrix[i,])
    }
    colnames(belief)[number_of_states+1] <- "belief_state"
    
    belief_proportions <- alpha-alpha
    colnames(belief_proportions) <- colnames(belief)[1:number_of_states]
    
    for (i in 1:dim(alpha_matrix)[1]) {
      c <- 0
      for (j in 1: dim(belief_matrix)[1]) {
        if (belief[j,ncol(belief)]==i) {
          c <-  c + 1
          belief_proportions[i,] <- belief_proportions[i,]+belief_matrix[j,]
        }
      }
      belief_proportions[i,] <- belief_proportions[i,]/c
    }
  } else {
    belief <- NULL
    belief_proportions <- NULL
  }
  
  ## ============Construct PG=====================
  if (is.null(belief_proportions)){
    pg <- NULL
  }else{
    belief_prop_matrix <- as.matrix(belief_proportions)
    pg <- create_PG(model = model,
          alpha = alpha,
          action_index_alpha = action_index_alpha,
          belief_matrix = belief_prop_matrix)
  }
  
  
  
  ### outputs and results
  
  ## producing the starting belief vector
  if (!is.character(start)) {
    if (sum(start)==1) {
      start_belief <- start
    }
  }
  # if the starting beliefs are given by a uniform distribution over all states
  if (sum(start== "uniform")==1) {
    start_belief <- rep(1/number_of_states,number_of_states)
  } else if (start[1]!="-") {  # if the starting beliefs include a specific subset of states
    # if the starting beliefs are given by a uniform distribution over a subset of states (using their names)
    if (!is.na(sum(match(start , states)))) {
      start_belief <- rep(0,number_of_states)
      start_belief[match(start,states)] <- 1/length(start)
    }
    # if the starting beliefs are given by a uniform distribution over a subset of states (using their numbers)
    if (!is.character(start)) { 
      if (sum(start)>=1 & length(start)<number_of_states) {
        start_belief <- rep(0,number_of_states)
        start_belief[start] <- 1/length(start)
      }
    }
  } else if (start[1]=="-") { # if the starting beliefs exclude a specific subset of states
    start_belief <- rep(1/(number_of_states-length(start)+1),number_of_states)
    if (is.na(as.numeric(start[2]))) {
      start_belief[match(start,states)] <- 0
    }
    if (!is.na(as.numeric(start[2]))) {
      start_belief[start] <- 0
    }
  }
  
  ## calculating the total expected reward for start_belief
  initial_belief_state <- which.max(alpha_matrix %*% start_belief)
  total_expected_reward <- max(alpha_matrix %*% start_belief)
  
  solution <- structure(list(
    method = method, 
    parameter = parameter,
    alpha = alpha,
    action_index_alpha = action_index_alpha,
    pg = pg,
    belief = belief, 
    belief_proportions = belief_proportions,
    total_expected_reward = total_expected_reward,
    initial_belief_state = initial_belief_state
  ), class = "POMDP_solution")
  
  structure(list(model = model,
    solution = solution,
    solver_output = solver_output
  ), class = "POMDP")
}

print.POMDP <- function(x, ...) {
  cat(
    "Solved POMDP model:", x$model$name, "\n", 
    "\tmethod:", x$solution$method, "\n",
    "\tbelief states:", nrow(x$solution$pg), "\n",
    paste0("\ttotal expected ", x$model$values, ":"), 
      x$solution$total_expected_reward,"\n\n" 
  )
}

print.POMDP_solution <- function(x, ...) {
 cat("POMDP solution\n\n")
 print(unclass(x))
}

5 Simulation of POMDP functions

5.1 Action_map function

# Action_map is used to get the action index for belief under policy
Action_map <- function(belief, policy = "optimal", alpha, action_index_alpha, obs){
  # if (policy == "fix"){
  #   if (belief[3] > 0.5 | belief[4] > 0.5){
  #     action_index <- 2
  #   }else if(belief[1] > 0.5 | belief[2] > 0.5){
  #     action_index <- 1
  #   }else{
  #     action_index <- 5 # action5==ignore
  #   }
  # }
  
  if (policy == "fix"){
    state_index <- which.max(belief)
    if (state_index == 3 | state_index == 4){
      action_index <- 2
    }else{
      action_index <- 1
    }
  }
  
  # if (policy == "fix"){
  #   if (belief[3] == 1 | belief[4] == 1){
  #     action_index <- 2
  #   }else{
  #     action_index <- 1
  #   }
  # }
  
  # =====fix policy based on observations====
  if (policy == "OSOM"){
    if (obs == 4 | obs == 5){
      action_index <- 2
    }else if(obs == 2 | obs == 3){
      action_index <- 1
    }else{ #if obs=1, observe nothing, do action5==ignore
      action_index <- 5
    }
  }

  # if (policy == "fix_1111"){
  #   if (obs == 1){
  #     action_index <- 5
  #   }else{
  #     action_index <- 1
  #   }
  # }
  
  if (policy == "optimal"){
    alpha_matrix <- as.matrix(alpha)
    get_index <- which.max(alpha_matrix %*% belief)
    action_index <- action_index_alpha[get_index]
  }
  
  if (policy == "random"){
    action_index <- sample(1:4, size = 1, replace = TRUE)
  }
  
  if (policy == "random_12"){
    action_index <- sample(1:2, size = 1, replace = TRUE)
  }
  
  if (policy == "e_greedy"){
    alpha_matrix <- as.matrix(alpha)
    get_index <- which.max(alpha_matrix %*% belief)
    action_index_opt <- action_index_alpha[get_index] #optimal action index
    
    action <- c(1:4)
    action_exclude <- action[action!= action_index_opt] #remove the optimal action index
    action_index_greedy <- sample(action_exclude, size = 1, replace = TRUE)
    
    action_index <- sample(c(action_index_opt, action_index_greedy),
                           size = 1, 
                           replace = TRUE,
                           prob = c(0.95, 0.05))
  }
  
return(action_index)
}

# Action_map(belief = c(0,0,0,1), 
#            policy = "fix",
#            alpha,
#            action_index_alpha)

5.2 Simulate_policy function

The main idea in this Simulate_policy function is to randomly generate \(state, observation\) based on their probabilities under \(action, belief\) at each time step.

Simulate_policy <- function(model_solved, n_trial = 1000, n_horizon = 60,
                                policy = "fix", b_0 = "random", seed = 123){
  
  states <- model_solved$model$states
  actions <- model_solved$model$actions
  observations <- model_solved$model$observations
  gamma <- model_solved$model$discount
  
  n_state <- length(states)
  n_obs <- length(observations)
  
  # policy
  if (policy == "optimal"){
    alpha <- model_solved$solution$alpha
    action_index_alpha <- model_solved$solution$action_index_alpha
  }
  
  if (policy == "e_greedy"){
    alpha <- model_solved$solution$alpha
    action_index_alpha <- model_solved$solution$action_index_alpha
  }
  
  # reward matrix
  reward <- model_solved$model$reward
  reward <- reward %>%
    mutate(start.state = fct_relevel(start.state, states))%>%
    spread(key = "action", value = "reward")%>%
    select(actions)
  
  # observations matrix
  observation_prob <- model_solved$model$observation_prob
  
  # transition matrix 
  transition_prob <- model_solved$model$transition_prob
  
  # # belief matrix
  # belief_matrix <- model_solved$solution$belief[, 1:n_state]%>%as.matrix()%>%unname()
  
  # if policy is fixed, add action5 = ignore
  if(policy == "fix" | policy == "OSOM"){
     #transition_prob[["ignore"]] <- diag(x= 1, nrow = n_state, ncol = n_state)
     # transition_prob[["ignore"]] <- matrix(data = c(1, 0, 0, 0,
     #                                                0, 1/3, 1/3, 1/3,
     #                                                0, 0, 0.5, 0.5,
     #                                                0, 0, 0, 1),
     #                                       nrow = n_state, ncol = n_state, byrow = TRUE)
     transition_prob[["ignore"]] <- matrix(data = 1/n_state,
                                           nrow = n_state, ncol = n_state, byrow = TRUE)
     # observation_prob[["ignore"]] <- matrix(data = c(1,0,0,0,0,
     #                                                 1,0,0,0,0,
     #                                                 1,0,0,0,0,
     #                                                 1,0,0,0,0),
     #                                        nrow = n_state, ncol = n_obs, byrow = T)

     observation_prob[["ignore"]] <- model_solved$model$observation_prob[[1]]
     reward <- reward%>%dplyr::mutate(ignore = .[,1]) #action5=ignore cost == action 1 cost

     model_solved$model$observation_prob <- observation_prob
     model_solved$model$transition_prob <- transition_prob
  }
  # 
  # ==========Simulation===============
  # run n_trial to get expected rewards
  # store results
  r_total_trials <- rep(0, n_trial)
  action_trials <- matrix(0, nrow = n_trial, ncol = n_horizon)
  state_trials <- matrix(0, nrow = n_trial, ncol = n_horizon)
  obs_trials <- matrix(0, nrow = n_trial, ncol = n_horizon)
  r_trials <- matrix(0, nrow = n_trial, ncol = n_horizon)

  set.seed(seed) # for reproduce
  #======outer loop: trials start=======
  for (t in 1:n_trial){
    discount <- 1
    r_total <- 0
    
    #====initial version 1: 
    if (b_0 == "uniform"){
      obs <- 1
      b <- rep(1/n_state, n_state)
      s_start <- sample(c(1:n_state), size = 1, replace = TRUE, prob = b) #random generate state: s_start
      }
    #===initial version 2: generate b_0 based on random observation at begining
    if (b_0 == "random"){
      obs <- sample(c(1:n_obs), size = 1, replace = TRUE)
      if (obs == 1){ #if observe nothing, randomly select b_0 from belief matrix
          # b_index <- sample(c(1:nrow(belief_matrix)), size = 1, replace = TRUE)
          # b <- belief_matrix[b_index, ]
          b <- rep(1/n_state, n_state)
        }else if(obs == 2){
          b <- c(1, 0, 0, 0)
          }else if(obs == 3){
            b <- c(0, 1, 0, 0)
            }else if (obs == 4){
              b <- c(0, 0, 1, 0)
              }else{
                b <- c(0, 0, 0, 1)
                }
      s_start <- sample(c(1:n_state), size = 1, replace = TRUE, prob = b) #random generate state: s_start
      }
    
    #====inner loop: periods start======
    for (h in 1:n_horizon){
      # get action index: a
      a <- Action_map(belief = b, policy = policy, alpha, action_index_alpha, obs)
      action_trials[t, h] <- a
      state_trials[t, h] <- s_start
      
      # random generate next state: s_end
      s_end <- sample(c(1:n_state), size = 1, replace = TRUE, prob = transition_prob[[a]][s_start, ])
    
      # random generate observation: obs
      obs_trials[t, h] <- obs
      obs <- sample(c(1:n_obs), size = 1, replace = TRUE, prob = observation_prob[[a]][s_end, ])
      
      # get the reward
      r <- reward[s_start,a]
      r_trials[t,h] <- r 
      r_total <- r_total + discount*r
      
      # update state, belief and discount 
      s_start <- s_end
      b_update <- belief_update(b, model = model_solved$model, act_index = a, obs_index = obs)
      b <- b_update
      discount <- discount*gamma
    }
    
    r_total_trials[t] <- r_total
  }
  
  return(list(r_total_trials = r_total_trials,
              r_trials = r_trials,
              action_trials = action_trials,
              state_trials = state_trials,
              obs_trials = obs_trials))
}

# Simulate_policy(model_solved = model_POMDPs_solved[[4]], 
#                                     policy = "fix",
#                                     n_trial = 1000,
#                                     n_horizon = 60)

5.4 Simulate_policy2 function

The main idea in this Simulate_policy2 function is to randomly generate \(state, observation\) at each time step, then the reward of each time step is \(\sum_s R(s,a)*belief(s)\).

# Simulate_policy2 <- function(model_solved, n_trial = 1000, n_horizon = 60,
#                                 policy = "fix", b_0 = "uniform"){
#   
#   states <- model_solved$model$states
#   actions <- model_solved$model$actions
#   observations <- model_solved$model$observations
#   gamma <- model_solved$model$discount
#   
#   n_state <- length(states)
#   n_obs <- length(observations)
#   
#   # policy
#   if (policy == "optimal"){
#     alpha <- model_solved$solution$alpha
#     action_index_alpha <- model_solved$solution$action_index_alpha
#   }
#   
#   # initial belief
#   if(b_0 == "uniform"){
#     b_0 <- rep(1/n_state, n_state)
#   }
#   
#   # reward matrix
#   reward <- model_solved$model$reward
#   reward <- reward %>%
#     mutate(start.state = fct_relevel(start.state, states))%>%
#     spread(key = "action", value = "reward")%>%
#     select(actions)
#   
#   # observations matrix
#   observation_prob <- model_solved$model$observation_prob
#   
#   # transition matrix 
#   transition_prob <- model_solved$model$transition_prob
#   
#   # ==========Simulation===============
#   # run n_trial to get expected rewards
#   r_total_trials <- rep(0, n_trial)
#   for (t in 1:n_trial){
#     discount <- 1
#     r_total <- 0
#     b <- b_0
#     s_start <- sample(c(1:n_state), size = 1, replace = TRUE, prob = b_0) #random generate state: s_start
#     for (h in 1:n_horizon){
#       # get action index: a
#       a <- Action_map(belief = b, policy = policy, alpha, action_index_alpha) 
#     
#       # random generate next state: s_end
#       s_end <- sample(c(1:n_state), size = 1, replace = TRUE, prob = transition_prob[[a]][s_start, ])
#     
#       # random generate observation: obs
#       obs <- sample(c(1:n_obs), size = 1, replace = TRUE, prob = observation_prob[[a]][s_end, ])
#     
#       # get the reward
#       r <- b %*% reward[,a]
#       r_total <- r_total + discount*r
#       
#       # update state, belief and discount 
#       b_update <- belief_update(b, model = model_solved$model, act_index = a, obs_index = obs)
#       b <- b_update
#       discount <- discount*gamma
#     }
#     
#     r_total_trials[t] <- r_total
#   }
#   return(r_total_trials)
# }

6 ASMBI: Multi-POMDPs

6.1 Step 1: Define and Run POMDPs for each type of hosts.

The defined POMDPs are stored in model_POMDPs and the solved results are stored in model_POMDPs_solved.

Run POMDPs: belief points are using “search” method to generate.

Run POMDPs: incremental pruning method

# # Define start belief, horizons, method, points type
# start<-"uniform"
# horizon <- NULL
# method <-"incprune" #method could be enum, twopass, linsup, witness, incprune, grid, mcgs
# params <-list(epsilon = 1e-1)  #set the stop tollerence 
# 
# # Construct POMDPs and Solve them
# ## initilize 
# model_POMDPs_2 <- rep(list(NULL), n_type_Hosts*n_error)
# model_POMDPs_solved_2 <- rep(list(NULL), n_type_Hosts*n_error)
# model_index <- 0
# 
# ## loop
# pb = txtProgressBar(min = 1, max = (n_type_Hosts*n_error), style = 3) #process bar 
# 
# for (err in 1:n_error){
#   for (o in 1:nOS){
#     for (c in 1:nCrit){
#       for (a in 1:nAdmin){
#       
#         #obtain the scenario of errors
#         Err_name <- errors_name[err] 
#         #obtain hosts' features
#         OS <- names(reward_list)[o]
#         Crit <- names(reward_list[[o]])[c]
#         Admin <- names(reward_list[[o]][[c]])[a]
#       
#         #obtain the corresponding transitions, observations and rewards
#         trans_prob <- transition_prob[[OS]]
#         obs_prob <- observation_prob[[OS]][[Crit]]
#         reward <- reward_list[[OS]][[Crit]][[Admin]]
#       
#         #consider errors in observation matrix
#         obs_prob <- Consider_Error_Observation(obs_matrix = obs_prob,
#                                              error_scan = errors_scan[err])
#       
#         #construct POMDP
#         model <- pomdp::POMDP(
#           name = paste(OS, Crit, Admin, Err_name,sep = "_"),
#           discount = discount,
#           states = states,
#           actions = actions,
#           observations = observations,
#           start = start,
#           transition_prob = trans_prob,
#           observation_prob =obs_prob,
#           reward = reward)
#         model_index <- model_index + 1
#         model_POMDPs[[model_index]]<-model
#         #progress bar
#         Sys.sleep(0.02)
#         setTxtProgressBar(pb, model_index)
#         #solve POMDP
#         model_solved_2 <- solve_POMDP2(model = model,
#                                     horizon = horizon,
#                                     method = method,
#                                     verbose = FALSE,
#                                     parameter = params)
#         model_POMDPs_solved_2[[model_index]] <- model_solved
#       }
#     }
#   }
# }

6.2 Step 2: Shape the outputs

  • Get the PG when considering 5 scenarios:
## Warnings in  Other_crit_FALSE_admin_FALSE_error_0 : The belief update is NAs because it is impossible to observe low observation under action auto-patching when the start belief is ( 0,0,1,0 )Warnings in  Other_crit_FALSE_admin_TRUE_error_0 : The belief update is NAs because it is impossible to observe low observation under action auto-patching when the start belief is ( 0,0,1,0 )Warnings in  Windows_crit_FALSE_admin_FALSE_error_0 : The belief update is NAs because it is impossible to observe low observation under action auto-patching when the start belief is ( 0,0,1,0 )Warnings in  Windows_crit_FALSE_admin_TRUE_error_0 : The belief update is NAs because it is impossible to observe low observation under action auto-patching when the start belief is ( 0,0,1,0 )Warnings in  Windows_crit_FALSE_admin_TRUE_error_0 : The belief update is NAs because it is impossible to observe low observation under action auto-patching when the start belief is ( 0,0,0,1 )Warnings in  Windows_crit_FALSE_admin_TRUE_error_0 : The belief update is NAs because it is impossible to observe medium observation under action auto-patching when the start belief is ( 0,0,0,1 )Warnings in  Windows_crit_TRUE_admin_TRUE_error_0 : The belief update is NAs because it is impossible to observe low observation under action auto-patching when the start belief is ( 0,0,1,0 )
  • Get the alpha vectors
  • Get the plot of policy graph

7 ASMBI-Simulation: Fixed Policy vs Optimal Policy

if observe nothing, taking action 5(ignore); if observe low and medium, taking action 1; if observe high and critical, taking action 2.

Here, we perform simulations for optimal, fix (1122), OSOM, random and epsilon greedy policeies. The number of horizons is 60. For each policy, the number of simulations is 1000. (The simulations are running for 60 horizons repeated 1000 times)

7.1 Simulation

## 
  |                                                                       
  |                                                                 |   0%
  |                                                                       
  |                                                                 |   1%
  |                                                                       
  |=                                                                |   1%
  |                                                                       
  |=                                                                |   2%
  |                                                                       
  |==                                                               |   3%
  |                                                                       
  |==                                                               |   4%
  |                                                                       
  |===                                                              |   4%
  |                                                                       
  |===                                                              |   5%
  |                                                                       
  |====                                                             |   6%
  |                                                                       
  |====                                                             |   7%
  |                                                                       
  |=====                                                            |   7%
  |                                                                       
  |=====                                                            |   8%
  |                                                                       
  |======                                                           |   9%
  |                                                                       
  |======                                                           |  10%
  |                                                                       
  |=======                                                          |  10%
  |                                                                       
  |=======                                                          |  11%
  |                                                                       
  |========                                                         |  12%
  |                                                                       
  |========                                                         |  13%
  |                                                                       
  |=========                                                        |  13%
  |                                                                       
  |=========                                                        |  14%
  |                                                                       
  |==========                                                       |  15%
  |                                                                       
  |==========                                                       |  16%
  |                                                                       
  |===========                                                      |  16%
  |                                                                       
  |===========                                                      |  17%
  |                                                                       
  |===========                                                      |  18%
  |                                                                       
  |============                                                     |  18%
  |                                                                       
  |============                                                     |  19%
  |                                                                       
  |=============                                                    |  19%
  |                                                                       
  |=============                                                    |  20%
  |                                                                       
  |=============                                                    |  21%
  |                                                                       
  |==============                                                   |  21%
  |                                                                       
  |==============                                                   |  22%
  |                                                                       
  |===============                                                  |  23%
  |                                                                       
  |===============                                                  |  24%
  |                                                                       
  |================                                                 |  24%
  |                                                                       
  |================                                                 |  25%
  |                                                                       
  |=================                                                |  26%
  |                                                                       
  |=================                                                |  27%
  |                                                                       
  |==================                                               |  27%
  |                                                                       
  |==================                                               |  28%
  |                                                                       
  |===================                                              |  29%
  |                                                                       
  |===================                                              |  30%
  |                                                                       
  |====================                                             |  30%
  |                                                                       
  |====================                                             |  31%
  |                                                                       
  |=====================                                            |  32%
  |                                                                       
  |=====================                                            |  33%
  |                                                                       
  |======================                                           |  33%
  |                                                                       
  |======================                                           |  34%
  |                                                                       
  |======================                                           |  35%
  |                                                                       
  |=======================                                          |  35%
  |                                                                       
  |=======================                                          |  36%
  |                                                                       
  |========================                                         |  36%
  |                                                                       
  |========================                                         |  37%
  |                                                                       
  |========================                                         |  38%
  |                                                                       
  |=========================                                        |  38%
  |                                                                       
  |=========================                                        |  39%
  |                                                                       
  |==========================                                       |  39%
  |                                                                       
  |==========================                                       |  40%
  |                                                                       
  |==========================                                       |  41%
  |                                                                       
  |===========================                                      |  41%
  |                                                                       
  |===========================                                      |  42%
  |                                                                       
  |============================                                     |  42%
  |                                                                       
  |============================                                     |  43%
  |                                                                       
  |============================                                     |  44%
  |                                                                       
  |=============================                                    |  44%
  |                                                                       
  |=============================                                    |  45%
  |                                                                       
  |==============================                                   |  45%
  |                                                                       
  |==============================                                   |  46%
  |                                                                       
  |==============================                                   |  47%
  |                                                                       
  |===============================                                  |  47%
  |                                                                       
  |===============================                                  |  48%
  |                                                                       
  |================================                                 |  48%
  |                                                                       
  |================================                                 |  49%
  |                                                                       
  |================================                                 |  50%
  |                                                                       
  |=================================                                |  50%
  |                                                                       
  |=================================                                |  51%
  |                                                                       
  |=================================                                |  52%
  |                                                                       
  |==================================                               |  52%
  |                                                                       
  |==================================                               |  53%
  |                                                                       
  |===================================                              |  53%
  |                                                                       
  |===================================                              |  54%
  |                                                                       
  |===================================                              |  55%
  |                                                                       
  |====================================                             |  55%
  |                                                                       
  |====================================                             |  56%
  |                                                                       
  |=====================================                            |  56%
  |                                                                       
  |=====================================                            |  57%
  |                                                                       
  |=====================================                            |  58%
  |                                                                       
  |======================================                           |  58%
  |                                                                       
  |======================================                           |  59%
  |                                                                       
  |=======================================                          |  59%
  |                                                                       
  |=======================================                          |  60%
  |                                                                       
  |=======================================                          |  61%
  |                                                                       
  |========================================                         |  61%
  |                                                                       
  |========================================                         |  62%
  |                                                                       
  |=========================================                        |  62%
  |                                                                       
  |=========================================                        |  63%
  |                                                                       
  |=========================================                        |  64%
  |                                                                       
  |==========================================                       |  64%
  |                                                                       
  |==========================================                       |  65%
  |                                                                       
  |===========================================                      |  65%
  |                                                                       
  |===========================================                      |  66%
  |                                                                       
  |===========================================                      |  67%
  |                                                                       
  |============================================                     |  67%
  |                                                                       
  |============================================                     |  68%
  |                                                                       
  |=============================================                    |  69%
  |                                                                       
  |=============================================                    |  70%
  |                                                                       
  |==============================================                   |  70%
  |                                                                       
  |==============================================                   |  71%
  |                                                                       
  |===============================================                  |  72%
  |                                                                       
  |===============================================                  |  73%
  |                                                                       
  |================================================                 |  73%
  |                                                                       
  |================================================                 |  74%
  |                                                                       
  |=================================================                |  75%
  |                                                                       
  |=================================================                |  76%
  |                                                                       
  |==================================================               |  76%
  |                                                                       
  |==================================================               |  77%
  |                                                                       
  |===================================================              |  78%
  |                                                                       
  |===================================================              |  79%
  |                                                                       
  |====================================================             |  79%
  |                                                                       
  |====================================================             |  80%
  |                                                                       
  |====================================================             |  81%
  |                                                                       
  |=====================================================            |  81%
  |                                                                       
  |=====================================================            |  82%
  |                                                                       
  |======================================================           |  82%
  |                                                                       
  |======================================================           |  83%
  |                                                                       
  |======================================================           |  84%
  |                                                                       
  |=======================================================          |  84%
  |                                                                       
  |=======================================================          |  85%
  |                                                                       
  |========================================================         |  86%
  |                                                                       
  |========================================================         |  87%
  |                                                                       
  |=========================================================        |  87%
  |                                                                       
  |=========================================================        |  88%
  |                                                                       
  |==========================================================       |  89%
  |                                                                       
  |==========================================================       |  90%
  |                                                                       
  |===========================================================      |  90%
  |                                                                       
  |===========================================================      |  91%
  |                                                                       
  |============================================================     |  92%
  |                                                                       
  |============================================================     |  93%
  |                                                                       
  |=============================================================    |  93%
  |                                                                       
  |=============================================================    |  94%
  |                                                                       
  |==============================================================   |  95%
  |                                                                       
  |==============================================================   |  96%
  |                                                                       
  |===============================================================  |  96%
  |                                                                       
  |===============================================================  |  97%
  |                                                                       
  |================================================================ |  98%
  |                                                                       
  |================================================================ |  99%
  |                                                                       
  |=================================================================|  99%
  |                                                                       
  |=================================================================| 100%

7.2 Shape the results

1. dataframe for rewards

2. dataframe for average rewards

3. dataframe for evolutions of states and actions

7.3 Plots

7.3.1 Costs

  1. For specific hosts under optimal policy, plot the costs boxplot across observation errors

  1. When error = 0, plot the costs boxplot over different policy

Exclue some results from the above plot:

  1. Average costs plot: error = 0% and exclude random policy

  1. Average costs plot for only optimal and OSOM: line plot over errors

  1. Plot only optimal and OSOM: error = 0%, bar plot

  1. Plot savings proportion: optimal vs OSOM when error = 0%

  1. Plot savings proportion: optimal vs fix

7.3.2 Plot the evolution of states

  • Here, we create Plot_evolution function to plot a heat map of evolutions for specific policy and type of hosts.
  • Plot_evolution_all: plot evolutions for type of host
Plot_evolution_all <- function(data_evolution, target = "states", trial_limit = 1000, plot_ratio = 0.1, 
                           OS = "Linux", 
                           Criticality = TRUE, 
                           Admin = TRUE,
                           Error = 0,
                           remove_random = TRUE,
                           remove_fix = TRUE){
  os <- OS
  crit <- Criticality
  admin <- Admin
  err <- Error
  limit <- trial_limit
  
  # symbolize
  target <- as.name(target)
  # enquotes
  target <- enquo(target)
  # string
  target_name <- quo_name(target)
  
  data <- data_evolution%>%
    filter(OS == os & Criticality == crit & Admin == admin & Error == err)
  
  policy <- data$policy
  
  data_list <- lapply(policy, function(p){
    data_p <- data%>%filter(policy == p)
    data_matrix <- data_p[[target_name]][[1]]
    data_df <- setNames(reshape2::melt(data_matrix), c("trial", "horizons", target_name))
    data_df <- data_df%>%mutate(policy = p)
    return(data_df)
  })
  
  data_all <- data_list%>%map_df(`[`) #convert list to dataframe
  
  if(remove_random & remove_fix){
    plot <- data_all%>%
    filter(trial <= limit & policy != "random" & policy != "fix")%>%
    ggplot(aes(x = horizons, y = trial, fill = !!target) )+
    facet_grid(~policy)+
    geom_tile(colour="white",size=0.01)+
    scale_fill_gradient2(low = "white", high = "black")+
    scale_y_continuous(expand=c(0,0))+
    scale_x_continuous(expand = c(0,0),
                       breaks = seq(0, 60, by = 6))+
    coord_fixed(ratio = plot_ratio)+
    theme_bw()+
    theme(axis.text.x = element_text(size = 5),
          legend.position = "bottom")
  }else if (remove_random & !remove_fix){
    plot <- data_all%>%
    filter(trial <= limit & policy != "random")%>%
    ggplot(aes(x = horizons, y = trial, fill = !!target) )+
    facet_grid(~policy)+
    geom_tile(colour="white",size=0.01)+
    scale_fill_gradient2(low = "white", high = "black")+
    scale_y_continuous(expand=c(0,0))+
    scale_x_continuous(expand = c(0,0),
                       breaks = seq(0, 60, by = 6))+
    coord_fixed(ratio = plot_ratio)+
    theme_bw()+
    theme(axis.text.x = element_text(size = 8),
          legend.position = "bottom")
  }else if(!remove_random & remove_fix){
    plot <- data_all%>%
    filter(trial <= limit & policy != "fix")%>%
    ggplot(aes(x = horizons, y = trial, fill = !!target) )+
    facet_grid(~policy)+
    geom_tile(colour="white",size=0.01)+
    scale_fill_gradient2(low = "white", high = "black")+
    scale_y_continuous(expand=c(0,0))+
    scale_x_continuous(expand = c(0,0),
                       breaks = seq(0, 60, by = 6))+
    coord_fixed(ratio = plot_ratio)+
    theme_bw()+
    theme(axis.text.x = element_text(size = 8),
          legend.position = "bottom")
  }else{
    plot <- data_all%>%
    filter(trial <= limit)%>%
    ggplot(aes(x = horizons, y = trial, fill = !!target) )+
    facet_grid(~policy)+
    geom_tile(colour="white",size=0.01)+
    scale_fill_gradient2(low = "white", high = "black")+
    scale_y_continuous(expand=c(0,0))+
    scale_x_continuous(expand = c(0,0),
                       breaks = seq(0, 60, by = 6))+
    coord_fixed(ratio = plot_ratio)+
    theme_bw()+
    theme(axis.text.x = element_text(size = 8),
          legend.position = "bottom")
  }
  
    
  return(plot)
}

Example:

  1. Plot evolution of states for all policies: Linux, critical=true, admin=true

  1. Plot evolution of actions for all policies: Linux, critical=true, admin=true

  1. Plot evolution of states for all policies: Linux, critical=true, admin=False

  1. Plot evolution of states for all policies: Windows, critical=true, admin=true

  1. Plot evolution of actions for all policies: Windows, critical=true, admin=true

  1. Compare Windows and Linux when policy = optimal

7.4 Observations (example)

Look at the state of each host over 21 month

Basically, we’re looking for the proportion of disappeared time in the 21 scanned period for each unique host. Then, we intend to obtain the mean of the proportion of disappeared time for each type of hosts.

For each unique host or IP, generate the vector of occurrence over months. Then, the proportion of missing is the number of 0’s/the totoal number of months. After that, we could get a set of proportions of missing for each OS type of hosts.

Then, we could draw the box plot for each OS type of hosts

Then, we can calculate the average, median of disappeared proportion and corresponding standard deviation for each group shown above.

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## 0.00000 0.09524 0.42857 0.45308 0.85714 0.95238